#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <quadmath.h>
#include <float.h>
#include <immintrin.h>

static __float128 truncate(__float128 x, int keep_bits)
{
   int e;
   frexpq(x,&e);
   __float128 scale=1.0;
   __float128 half=0.5;
   int shift=0;
   if ( e-1 < keep_bits)
   {
      shift= keep_bits-(e-1);
      scale= ldexpq(half, shift+1);
   }
   if ( x>0.0)
      return floorq(x*scale)/scale;
   else
      return ceilq(x*scale)/scale;
}

static inline uint64_t as_uint64(double f)
{
  union
  {
    double f;
    uint64_t i;
  } u = {f};
  return u.i;
}

static void print_lines(const char* lines[], int line_count)
{
   for (int i=0;i<line_count;i++)
      printf("%s\n",lines[i]);
   printf("\n");
}

static void print_first_part()
{
   const char* lines[]={
      "#include <stdint.h>",
      "",
      "#define DBL_EXP_BASE 1023",
      "",
      "#define  likely(x)  __builtin_expect((x), 1)",
      "#define  unlikely(x)  __builtin_expect((x), 0)",
      "",
      "#define C2 (1.0/(1.0*2.0))",
      "#define C3 (1.0/(1.0*2.0*3.0))",
      "#define C4 (1.0/(1.0*2.0*3.0*4.0))",
      "#define C5 (1.0/(1.0*2.0*3.0*4.0*5.0))"
   };
   print_lines(lines, sizeof(lines)/sizeof(lines[0]));
}

static void print_external_varibles()
{
   const char* lines[]={
      "// exp_data_52[i]= (2^(i/128)-1.0)*2^52",
      "uint64_t exp_data_52[128];",
      "",
      "// exp_data_61[i]= 2^(i/128)*2^61",
      "uint64_t exp_data_61[128];",
      "",
      "// exp_data_64[i]= (2^(i/128)-1.0)*2^64",
      "uint64_t exp_data_64[128];"
   };
   print_lines(lines, sizeof(lines)/sizeof(lines[0]));
}

void print_2_double_const()
{
   char buf[64];   
   __float128 q2  = 2.0Q;
   __float128 tmp;
   __float128 ln2 =  M_LN2q;
   
   quadmath_snprintf (buf, sizeof buf, "%.20Qf", ln2);
   printf("#define LN2\t%s\n",buf);
   quadmath_snprintf (buf, sizeof buf, "%.20Qf", 1.0Q/ln2);
   printf("#define INV_LN2\t%s\n",buf);

   __float128 b  = ln2;
   __float128 b0 = truncate(b,36);
   printf("#define Ln2hiN  (%la)\n",(double)b0);
   printf("#define Ln2loN  (%la)\n\n",(double)(b-b0));

   int e;
   double f,eps;
   uint64_t u64_n;

   int  ns[2]={961,1022};
   for (int i=0;i<2;i++)
   {
      tmp= ln2 * (__float128)ns[i];
      f=(double)tmp;
      frexp(f,&e);
      eps = ldexp(0.5,e-52);
      u64_n= as_uint64(f-eps);
      printf("// LN2_X_%d_AS_UINT64=ln(2)*%d\n", ns[i],ns[i]);
      printf("#define LN2_X_%d_AS_UINT64   0x%016" PRIx64 "\n\n",ns[i], u64_n);
   }
}

int main(int argc, char* argv[])
{
   print_first_part();
   print_2_double_const();
   print_external_varibles();
   return 0;
}

